One-dimensional diffuse scattering of 1,3-di(tert-butyl)cyclopentadienyl pentaphosphaferrocene modelled with closed-form expressions

The diffuse scattering of 1,3-di(tert-butyl)cyclopentadienyl pentaphosphaferrocene (Cp′′FeP5) is modelled with closed-form expressions derived from a growth model with the range of interactions (reichweite) s = 2.


Introduction
Disordered modular structures can be seen as a generalized form of crystalline matter where the arrangement of the individual (ordered) modules (layers, rods, bricks) is ambiguous and governed by probabilities. The case of layer structures with stacking disorder is well understood, since it can be described with simple growth models. In such models, the orientation of a layer depends on the orientation of a finite number s ! 1 [the range of interaction (Treacy et al., 1991) or reichweite (Kakinoki & Komura, 1954)] of prior layers. In this context, the term 'growth' is to be understood formally with respect to the model, not the growth of a given crystal. The stacking arrangement of the crystal may also have been formed during a phase transition or by ageing.
The diffuse scattering caused by layer structures following a growth model has been intensely investigated [see e.g. Treacy et al. (1991) or Welberry (2010) and references therein]. Perhaps surprisingly, the function series used to calculate diffuse scattering are distinctly better behaved than the corresponding series of periodic structures, which typically do not converge at any point.
In particular, owing to pointwise exponential convergence, often only few layers have to be taken into account to adequately simulate diffuse scattering. In a sense, after these few layers, the 'average' structure, a weighted overlay of all the stacking possibilities, is obtained. This is the principle implemented in general programs such as DIFFaX (Treacy et al., 1991).
However, in many common cases, instead of infinite series also closed-form expressions can be derived [see e.g. Jagodzinski (1949) and Kakinoki & Komura (1954)]. On the one hand, these are less general, as they are limited to special cases. On the other hand, they might be computationally more favourable and show no convergence problems in the highly correlated case. Moreover, they can be analytically differentiated for use in least-squares optimization problems. Above all, closed-form expressions may provide more insight into the shape of the diffraction peaks or the reason for homometry.
Recently, we described crystallization experiments of the ferrocene analogue Cp 00 Fe( 5 -P 5 ), where one cyclopentadienyl (Cp, 5 -C 5 H 5 À ) ring has been formally replaced by Cp 00 = 5 -t Bu 2 C 5 H 3 À , the 1,3-di-tert-butyl substituted analogue of Cp and the second Cp ring has been formally replaced by an aromatic 5 -P 5 À ring. Different crystals extracted from the same sample featured different degrees of disorder (Peresypkina et al., 2022). One of the crystals featured especially pronounced diffuse scattering with reflections indicating two different polytypes. Here, the diffuse scattering will be explained qualitatively and modelled quantitatively with growth model-derived closed-form expressions.

Diffraction
Diffraction intensities from crystals of Cp 00 FeP 5 prepared according to Fleischmann et al. (2015) were collected at the P24 beamline of the PETRA III synchrotron, DESY (Germany) at 20 K (Peresypkina et al., 2022). For all crystals, the structure was refined in the Pca2 1 space group with more or less disorder of the Cp 00 FeP 5 molecule about the y = 1 4 pseudo-reflection plane. For the crystal described below, the occupation of the two orientations was refined to 72.3:27.7 (3). The minor component was ignored for the modelling of the diffuse scattering. For further data collection and refinement details as well as atomic coordinates, see Peresypkina et al. (2022). Atomic form factors where approximated using the polynomial coefficients tabulated in Brown et al. (2006).
Reciprocal space sections were calculated using the CrysAlisPro software (Rigaku Oxford Diffraction, 2021). One-dimensional profiles were then extracted by summing up over 20 pixels perpendicular to the rod. The centres of the rods were determined visually.

Stacking arrangements
The possible stacking arrangements have already been discussed in Peresypkina et al. (2022) and will be briefly recapitulated here. The structure is composed of ordered monoclinic/rectangular layers with p1a1 symmetry and one molecule located on the general position ( Figs. 1 and 2). The translation lattice (henceforth simply lattice) of the layers is spanned by (a, b). The width of the layer is defined by the length of the c 0 vector perpendicular to the layer plane.
Given the n th layer, the (n + 1) st layer can appear in one of two positions, which is generated from the n th layer by application of a glide reflections with the intrinsic translations c 0 or b/2 + c 0 . These two operations will be called the c-and nglide for convenience and are indicated using the corresponding graphical symbols in Fig. 2. The c-glide plane is located at x = 0 and the n-glide plane at x = 1 4 . According to these rules, the origin of the n th layer can appear either at nc 0 or at a/2 + b/2 + nc 0 . Thus, every polytype can be described by a family of integers ð n Þ n2Z with n = 0, 1. The origin of the n th layer then is n (a + b)/2 + nc 0 . If two adjacent layers n and n + 1 are related by a c-glide then n+1 = n , if they are related by a n-glide then n+1 = 1 À n .

MDO polytypes
In the context of OD theory, which describes polytypes that are locally equivalent, polytypes of a maximum degree of order (MDO) play a special role. The MDO polytypes of a family cannot be decomposed into fragments of simpler The Pca2 1 structure of Cp 00 FeP 5 viewed down [010]. Fe (dark orange), P (bright orange) and C (grey) atoms are represented by ellipsoids drawn at the 50% probability level, H atoms by white spheres of arbitrary radius. Dotted curves mark the interface between layers.

Figure 2
A single layer of Cp 00 FeP 5 molecules projected on the layer plane (001). The a-glide planes of the p1a1 layer group are indicated by the conventional symbols (Hahn & Aroyo, 2016) in black. Potential glide planes relating the layer to the next one are given in red. The symbols in this case are to be read with respect to the (a, b, 2c 0 ) basis. H atoms are omitted for clarity. polytypes, i.e. polytypes of only a subset of n-tuples of adjacent layers (Dornberger-Schiff, 1982). MDO polytypes can be considered as the 'alphabet' of an OD polytype family: all polytypes can be decomposed into fragments of MDO polytypes. Moreover, in the majority of cases, ordered polytypes belong to the MDO class.
Cp 00 FeP 5 forms non-OD polytypes because pairs of adjacent layers related by c-or n-glides are not equivalent. However, the MDO concept is just as useful in this case. Here, there are two MDO polytypes: Pca2 1 [all adjacent layers related by cglides, ( n ) = . . . , 0, 0, 0, 0, . . . ] and Pna2 1 [all adjacent layers related by n-glides, ( n ) = . . . , 0, 1, 0, 1, . . . ]. In both cases c = 2c 0 . The two MDO polytypes are shown in Fig. 3. Note that the 2 1 screw rotations are at different positions with respect to the layers. On rods h + k even only sharp reflections are observed, whereas on rods h + k odd, distinct diffuse scattering is apparent. The positions of all maxima are in agreement with a lattice spanned by (a, b, 2c 0 ), which corresponds to either of the two MDO polytypes. Qualitatively, the polytypes can be differentiated by the systematic absences of the c [100] or n [100] glide reflections, respectively. For 0kl reflections where k is odd (for k even, both glide reflections feature the same reflection conditions), strong reflections l odd suggest the Pca2 1 polytype. Additional weak reflections l even prove a non-negligible contribution of Pna2 1 fragments [ Fig. 4(a)].

Qualitative interpretation of the diffraction pattern
Thus, it appears that the crystal is built of both polytypes. An intergrowth of two or more distinct polytypes is called an allotwin (Nespolo et al., 1999). However, in contrast to the crystal described here, the domains of an allotwin are macroscopic.
Additional extremely weak reflections are observed at halfintegral k values, which are due to superstructure formation in the [010] direction for the Pca2 1 polytype as described in Peresypkina et al. (2022). These reflections are significantly more pronounced for pure Pca2 1 crystals. But even there the modulation is minute. For the disordered crystal described here, these reflections can be neglected without hesitation.

Quantitative interpretation of the diffraction pattern
In the following, positions in reciprocal space will be expressed with respect to the reciprocal basis ða Ã ; b Ã ; c Ã 0 Þ = ða=jaj 2 ; b=jbj 2 ; c 0 =jc 0 j 2 Þ. The coordinates are given as hk to emphasize that diffraction intensities can only appear at integral h and k owing to the layer lattices, but diffuse scattering may appear along rods parallel to c Ã 0 at arbitrary real values, in the case of disordered layer arrangements.
Let F n (hk), h; k 2 Z, 2 R be the structure factor of the n th layer. Layers n even are translationally equivalent and therefore their structure factor only differs in phase from F 0 (hk). In contrast, layers n odd are reflected at [100] and thus will be expressed with respect to F À 0 ðhkÞ ¼ F 0 ðhkÞ. For brevity, henceforth the argument hk of F 0 and F À 0 will be omitted. The structure factor of the n th layer is given by: The overall diffraction intensity of a structure can be written as where an asterisk indicates the complex conjugate and Án designates the distance between layers. By separating into even and odd Án and substituting equation (1): For h + k even, ( n À n+1 )(h + k)/2 is integral and therefore I(hk) is independent of ( n ). Thus, all polytypes produce the same intensity distribution on rods h + k even. In particular, by application of the Dirichlet kernel, one can show that intensities appear only for integral and half-integral : The term expð2iÞ is 1 and À1 for integral and half-integral , respectively. These reflections are called family reflections, since they are identical for all members of the polytype family.

The special case of the h = 0 plane
Qualitatively, the diffuse scattering is best understood for the h = 0 plane. As noted above, the Pca2 1 and Pna2 1 MDO polytypes can be immediately distinguished by reflections on the k odd rods. The former features sharp reflections at integral, the latter at half-integral values.
Quantitatively, the situation is likewise simpler, because F(0k) = F À (0k) and equation (9) simplifies to In the simplest growth model, the probability of a c-glide relating adjacent layers is given by P (and conversely the probability of an n-glide by 1 À P). However, such a simple model cannot explain the existence of both, the Pca2 1 and Pna2 1 reflections. In fact, these models produce on rods h + k odd streaks of the form where c = 2P À 1 and d c stands for the family of functions These functions are well known and generally describe the shape of nearest-neighbour (s = 1) models with the nearestneighbour correlation c (Welberry, 2010;Stö ger et al., 2021). For c approaching 1 and À1, d c () converges to a Dirac comb at integral and half-integral , respectively. With decreasing |c|, the peaks become successively broader and at c = 0, d c () is constant. However, according to this model, in no case peaks are observed simultaneously for integral and half-integral , which is in contradiction with the observed intensities (Fig. 5).
Thus, a range of interaction s ! 2 is required. In an s = 2 two-neighbour growth model, the state of the previous step is considered as in the Markov chain: The parameter A gives the probability that a c-glide follows a c-glide and B that an n-glide follows an n-glide. For example, A = 1, B < 1 and A < 1, B = 1 represent the Pca2 1 and Pna2 1 polytypes, respectively and A = B = 0 is the ordered non-MDO polytype . . . cncn . . . . For A = 1 À B, the model degenerates to the single-neighbour model P = A = 1 À B described above.
Henceforth, degenerate Markov chains A = 1 or B = 1 as well as the ordered A = B = 0 model will be disregarded. Then, it is easy to show that the chain converges to an equilibrium state where the probability of c-and n-glides is and To calculate the diffuse scattering of a given model, the probabilities P Án ¼ Prð n ¼ nþÁn Þ have to be calculated [see equation (9)]. Since n can adopt two values, the number of states is doubled to four: The state c0 means a c-glide after n = 0, etc. The n th state of the Markov chain will be described by the row-vector n ¼ ½Prðc0Þ; Prðn0Þ; Prðn1Þ; Prðc1Þ: According to equations (14) and (15), the initial state of the chain is The transition probability matrix P, which gives the progression of the chain according to nþ1 ¼ n P is Therefrom, the diffuse scattering on rods h + k odd can be calculated as where For derivation see Appendices A1 and A2 and also Kakinoki & Komura (1954). Two cases can be distinguished. For large A or B (white region in Fig. 6 Sign of Q depending on A and B: white '+', blue 'À' and black 0.  correlations and the weighting of the two terms is not immediately obvious, a qualitative analysis confirms the expectations: for large A and small B, sharp diffraction spots are observed at integral , corresponding to the Pca2 1 MDO polytype [red curve in Fig. 7(a)]. With increasing B, additional reflections slowly appear at half-integral , corresponding to the Pna2 1 polytype (black and blue curve). Ultimately, for A = B the diffraction pattern corresponds to that of a mix of Pca2 1 and Pna2 1 polytypes with equal volume weights and peaks enlarged along c*.
With small A and B (blue region in Fig. 6), Q is negative and therefore (Q) 1/2 and R are purely imaginary. The case is shown in Fig. 7(b) for small A. For large B 'reflections' of the Pna2 1 polytype are observed (green curve). With decreasing B, these reflections split (black and blue curve). Ultimately, when B approaches A (red curve), reflections are observed at uneven fourths of , which corresponds to the . . . cncn . . . non-MDO polytype with c = 4c 0 , leading to the doubling of the cell parameter c with respect to the experimentally observed polytypes. Reflections at integral and half-integral values are not observed, because this particular polytype features Pna2 1 symmetry.
The problem of the imaginary value of (Q) 1/2 in the case of Q < 0 can by solved by a purely real form of equation (23), given in Appendix A2. It might be more appropriate for calculations, but since the diffraction pattern from the crystal does not suggest this case, it is not discussed any further.
For Q approaching 0 (black line in Fig. 6), the denominator of R goes to 0 and R diverges. However, in that case c 1 = c 2 = c = (A À B)/2 and R can be set to an arbitrary value (such as R = 0) to get IðhkÞ / jF 0 ðhkÞj 2 d c ðÞ ð 25Þ which is the diffraction pattern of a nearest-neighbour model with correlation c = (A À B)/2 [see equation (11)].
This kind of homometry is well known (Welberry, 2010) and has to be taken into account when determining the parameters of the growth model. In fact, homometry or quasi-homometry may lead to local minima and therefore a local search may not be sufficient.   preference for the former. Thus, quantitatively one would expect probabilities A > B > 1 2 . To derive A and B, a rod h = 0, k odd was chosen where the diffraction intensities were strong and the peak maxima were consistent with the integrated diffraction intensities. The (01)* rod was the optimal rod in that regard.

Estimation of the parameters
The high-quality diffraction data obtained with the highflux non-divergent synchrotron beam and a noiseless direct photon-counting detector at the P24 beamline (DESY, Hamburg) operated in shutterless mode allowed very narrow scans (0.1 2). Thus, little experimental broadening is expected, which makes it possible to reliably reconstruct and quantify the diffuse scattering.
The experimental peak broadening was approximated by a Gaussian distribution, whose variance 2 was determined by least squares refinement against the sharp family reflections (02l) (Fig. 9).
Intensities of the (01)* rod were calculated using the analytical expression equation (20) and numerically convoluted with the Gaussian peak broadening. The scale factor was determined by linear regression after each calculation.
The origin ( = 0) and the length |c Ã 0 | in pixels as well as A and B were determined using a combination of global search with local searches according to the multi-coordinate search (MCS) approach (Huyer & Neunmaier, 1999) (Fig. 10) implemented in custom routines. The minimized function was the square of the difference between measured and calculated intensities with unit weight: Since each rod is calculated in ms times on a standard computing hardware, the search ended in a short time with a final residual of R p = 1.2 %. Only one 'basket entry' (Huyer & Neunmaier, 1999) was obtained, which means that the function had only a single local minimum. Thus, homometry was not a problem in this case.
The refined values are summarized in Table 1, left row. As predicted, c-glides induce c-glides (83.1% probability) and nglides induce n-glides (61.1% probability). In that sense, as we had expected, the crystal can be considered as a disordered equivalent of an allotwin. In contrast to the crystal under investigation, though, an allotwin is composed of macroscopic domains, i.e. A and B are very close to 1.
Since the c-glide induces the n-glide with lower probability than the n-glide induces the c-glide (16.9% versus 38.9%), the Pca2 1 fragments prevail in agreement with the experimentally observed systematic absences and relative weight of the disordered components in structure refinements.
3.8. The general case of planes h 6 ¼ 0 The h 6 ¼ 0 planes are more difficult to interpret quantitatively, because the characteristic reflections produced by both MDO polytypes overlap. The simplification F(0k) = F À (0k) does not apply and the general equation (9) has to be used. The intensity distribution for rods h + k odd then calculates as Measured intensities of the (02)* rod and calculated profile using Gauss functions with = 2.96646 pixels.

Figure 10
Measured and simulated intensities of the (01)* rod.
For derivation see Appendix A3.
The shape functions of the first term are weighted by the average of the auto-correlations F 0 F Ã 0 and F À 0 F ÀÃ 0 . Note that here, the correlations c 1 and c 2 enter as their squares and the argument of the shape function d is 2 instead of . This is reasonable, because the term considers only layer pairs distanced by an even number Án of layers.
The shape functions of the second term are weighted by the average cross-correlations F 0 F ÀÃ 0 and F À 0 F Ã 0 . This average is real, but may be negative. However, the sum of auto-and cross-correlations is , which ultimately prevents negative intensities. Owing to the two different terms in equation (27), the peaks can adopt a distinctly skewed shape.
To exemplify equation (27), Fig. 11(a) gives the sum of autocorrelations F 0 F Ã 0 þ F À 0 F ÀÃ 0 , the sum of the cross-correlations F 0 F ÀÃ 0 þ F À 0 F Ã 0 and the intensity of a (21)* rod for the parameters A = 0.85 and B = 0.75. Fig. 11(b) plots the profiles of the reflections of equation (27) under the assumption of constant and equal F 0 and F À 0 and the dissection into the two contributing terms. The result of a MCS refinement using equation (27) is shown in Fig. 12 and summarized in Table 1. The fit is still reasonable, though slightly worse than in the (01)* case.

Significance of refinements
Comparing the refinements on the (01)* and (21)* rods, the A values agree well. However, a distinctly higher B value was derived from the (21)* rod. Indeed, since the contribution corresponding to the Pna2 1 MDO polytype is distinctly   (27) if F 0 and F À 0 were constant and equal. The sum of the two terms is also shown.

Figure 12
Measured and simulated intensities of the (21)* rod.   Table 1. R p values are encoded by colours on a logarithmic scale. less pronounced (R > 0), the B parameter is defined worse than the A parameter. This can be seen by plotting the loss function against the A and B parameters with fixed metric parameters. The minimum lies in a valley that is much steeper in A than in B direction.
The difference in the B value cannot however be explained by imprecision alone, since the valleys are located at distinctly different positions (Fig. 13), systematically shifted toward the larger B values for the (21)* rod.
The stability of A and the higher B for h 6 ¼ 0 is confirmed by additional refinements against the (03)* an (23)* rods (Table 2). These refinements are less reliable, because the intensities on these rods are significantly weaker and the contribution of B is even less significant.
For k = 5 rods, the extracted intensities are too inconsistent with the structure factors for reasonable refinements. Likewise for h = 1 rods with even k contribution of B to the shape was too small for a meaningful determination of this value.
One source of error might be inadequate intensity extraction. Indeed, the relative peak heights sometimes disagree with the integrated I hkl values. However, we suppose that the biggest error is due to neglecting desymmetrization. A layer can appear in three unique environments with respect to the two neighbours: cc, nn and cn (= nc). However, structural data only exists for the cc (corresponding to the Pca2 1 MDO polytype) case. For improved simulations, one would either have to grow crystals featuring the other layer contacts or resort to theoretical structure optimizations.
For the (01)* rod this error seems to be minimized. In fact, the FF À* + F*F À cross-correlation term does not appear in the analytical expressions, since in projection along [100] both orientations are identical. We suppose that also the desymmetrization is less pronounced in the [100] projection, leading to a more reliable estimation of B. Moreover, in that case, the values of A and B are derived from distinct 'reflections', whereas for h 6 ¼ 0, they are 'encoded' in the same reflections.

Conclusion and outlook
Having high-quality diffraction data for single crystals of Cp 00 FeP 5 obtained with high-flux synchrotron radiation at the P24 beamline (DESY, Hamburg), we attempted to describe the shape of the experimentally observed one-dimensional diffuse scattering with closed-form expression derived from a growth model with range of interaction s = 2 and explain the average disorder in the structure (Peresypkina et al., 2022). Thus, a direct relation between the stacking-fault probabilities of both MDO polytypes and the form of the diffraction maxima could be drawn.
Closed-form expressions allow for very fast calculations of diffuse scattering, which we used for respective refinements. The fact that tens of thousands of rods can be calculated in seconds enabled a global search. This is crucial in the case of diffuse scattering, since homometry may lead to multiple local minima. Moreover, closed-form expressions directly explain the shape and position of peaks in the diffuse scattering and allow a direct derivation of homometric pairs of disordered stacking arrangements.
However, such an approach is not general. Currently, deriving the expressions for diffuse scattering intensity is tedious and error-prone. In the future this might be automatized by symbolic algebra, since the theory of Markov chains is well understood. Even for the title compound, owing to the different orientations of the even and odd numbered layers, the general expressions are rather unwieldy. These complications disappeared, when only considering the (0k)* plane. Effectively, this means looking at a projection along [100], for which the orientation with respect to [100] vanishes. Then, the diffuse scattering can be described as the sum of two nearest-neighbour models. In the general case (h 6 ¼ 0), however, the expressions become less intuitive. With an increasing range of interaction or number of orientations, the mathematical expressions will become more and more cumbersome. A fundamental complexity limit is due to tha fact that, in general, only roots of polynomials up to degree four can be expressed by radicals. Thus, for more complex problems, only numerical solutions can be given.
The imperfect fit is most likely due to ignoring desymmetrization effects. We will show such a case in an upcoming publication of a different molecule, where two kinds of polytypes could be grown and thus desymmetrization could be taken into account. When considering different environments, Markov chains with yet more states have to be used.

A1. Correlations of the s = 2 model
Given an initial state 0 , the n th state n is calculated by n ¼ 0 P n . To calculate P n , the matrix P is diagonalized, i.e. expressed as QDQ À1 , where D is a diagonal matrix (Kakinoki & Komura, 1965). This reduces the problem to calculating P n ¼ QD n Q À1 , which is trivial, since the n th power of a diagonal matrix is obtained by taking the n th power of the diagonal elements. Matrix diagonalization is an eigenvector and eigenvalue problem. For P of equation (19) we obtain (as one possible solution) where Given the probabilities of c-and n-glides according to equations (14) and (15), the initial state of the chain is and n th state is However, we are not interested in the full state only the sum of the states for which = 0, 1, that is P Án ¼ n ð1Þ þ n ð2Þ. By substitution of equations [ (29) and (30)] into equation (33), one obtains: where Thus, in general convergence to the equilibrium state lim n!1 P Á2n ¼ 1 2 is given by the sum of two power series. The correlation c Án = 2P Án À 1 is: Note that if Q is negative, (Q) 1/2 and R are purely imaginary. However, in that case 1 + R and 1 À R are complex conjugate, as c 1 and c 2 are. Thus, c Án is the sum of two complex conjugates and therefore real, as expected.
A2. Diffuse scattering for h = 0 Substituting of equation (38) into equation (10, one obtains: The equivalence of the sums over Án and the shape functions d c 1 and d c 2 have for example been derived in Stö ger et al. . For negative Q, the terms in equation (40) are complex. Since both terms are complex conjugate, resulting in a positive real intensity, as expected. By using the identities ðc 1 Þ 2 þ ðc 2 Þ 2 ¼ ðA À BÞ 2 þ 2ðA þ B À 1Þ; ðc 1 Þ 2 c 2 þ c 1 ðc 2 Þ 2 ¼ ÀðA À BÞðA þ B À 1Þ: Equation (40) can also be written using only real terms as where C = A + B À 2. In the Q < 0 case, this can not generally be interpreted in terms of contributions by the Pca2 1 and Pna2 1 polytypes. One exception is A = B, in which case equation (45) becomes where c = 2A À 1 = 2B À 1.
A3. Diffuse scattering for h 6 ¼ 0 Substituting of equation (38) into equation (9), one obtains: Here, d 0 c is a family of shape functions that can be derived in analogy to d c :